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The magnetic susceptibility of the two-dimensional repulsive Hubbard model with nearest- 
neighbor hopping is investigated using the diagram technique developed for the case of strong 
correlations. In this technique a power series in the hopping constant is used. At half-filling the 
calculated zero- frequency susceptibility and the square of the site spin reproduce adequately results 
of Monte Carlo simulations. Also in agreement with numerical simulations no evidence of ferromag- 
netic correlations was found in the considered range of electron concentrations 0.8 ^ n ^ 1.2 for 
the repulsion parameters 8\t\ < U < 16\t\. However, for larger U/\t\ and |1 — n| « 0.2 the nearest 
neighbor correlations become ferromagnetic. For n <■ 0.94 and n ^ 1.06 the imaginary part of the 
real-frequency susceptibility becomes incommensurate for small frequencies. The incommensurabil- 
ity parameter grows with departure from half-filling and decreases with increasing the frequency. 
This behavior of the susceptibility can explain the observed low-frequency incommensurate response 
observed in normal-state lanthanum cuprates. 

PACS numbers: 71.10.Fd, 71.27.+a, 75.40.Gb 



I. INTRODUCTION 

The Hubbard modelA is thought to be appropriate to 
describe the main features of electron correlations in nar- 
row energy bands, leading to collective effects such as 
magnetism and metal-insulator transition. It has been 
often used to describe real materials exhibiting these phe- 
nomena (see, e.g., Refs.[3,[l, and references therein). 

In more than one dimension, the model is not exactly 
solvable and a variety of numerical and analytical ap- 
proximate methods was used for its study. Among oth- 
ers there are Monte Carlo simulationsj^^ different clus- 
ter methods^ the composite operator formalism)^ the 
generating functional approach,^ Green's function de- 
coupling schemes,^ and variational approaches. Along 
with these methods various versions of the diagram 
technique-i^iiiiii^ii^i^ have been used for the investiga- 
tion of the model. In the case of strong electron corre- 
lations when the ratio of the hopping constant t to the 
on-site repulsion [/ is a small parameter the use of the 
diagram technique based on the series expansion in this 
parameter is quite reasonable. 

In the present work we use the diagram technique of 
Refs. [12 and [ij for investigating the magnetic suscep- 
tibility of the one-band two-dimensional repulsive Hub- 
bard model with nearest-neighbor hopping in the case 
of strong electron correlations. In this version of the di- 
agram technique terms of the power expansion are ex- 
pressed through cumulants of creation and annihilation 
electron operators. The considered model possesses the 
electron-hole symmetry and results obtained for electron 
concentrations n < 1 are replicated for n> I. Therefore 
in the following discussion we shall restrict our consider- 
ation to the former region of concentrations. 



We found that at half-filling the calculated tempera- 
ture dependence of the zero-frequency susceptibility re- 
produces adequately key features of results of Monte 
Carlo simulationsi^ The uniform susceptibility tends to 
a finite value for vanishing temperature. The stag- 
gered susceptibility diverges with decreasing temperature 
which signals the establishment of the long-range anti- 
ferromagnetic order. The transition temperature Tq is fi- 
nite which indicates the violation of the Mermin- Wagner 
theoremji^ However, the transition temperature is al- 
ways lower than the analogous temperature in the ran- 
dom phase approximation (RPA). Besides, the transition 
temperature decreases with decreasing the ratio \t\/U of 
the hopping constant and the on-site repulsion, i.e. the 
violation of the Mermin- Wagner theorem becomes less 
pronounced on enforcing the condition for which the ap- 
proximation was developed. For small ratios \t\/U the 
calculated square of the site spin differs by less than 10% 
from the data of Monte Carlo simulations. Also in agree- 
ment with Monte Carlo results we found no evidence 
of ferromagnetic correlations in the considered range of 
electron concentrations 0.8 ^ n ^ 1.2 for the repul- 
sion parameters 8|t| < U < 16\t\. However, for larger 
U/\t\ and |1 — n| w 0.2 the nearest neighbor correla- 
tions become ferromagnetic. In the case U = 8\t\ for 
0.94 ^ n < 1.06 the zero-frequency susceptibility and 
the imaginary part of the susceptibility for low real fre- 
quencies are peaked at the antiferromagnetic wave vec- 
tor (tt, tt). For smaller and larger concentrations these 
susceptibilities become incommensurate - momenta of 
their maxima deviate from (tt, tt) - and the incommensu- 
rability parameter, i.e. the distance between (tt, tt) and 
the wave vector of the susceptibility maximum, grows 
with departure from half-filling. With increasing the fre- 




FIG. 1: The diagram equation for Z)(k, ioj,/). 



quency the incommensurability parameter decreases and 
finally vanishes. This behavior of the strongly correlated 
system resembles the incommensurate magnetic response 
observed in the normal-state lanthanum cuprates^^ and 
can be used for its explanation. 

Main formulas used in the calculations are given in the 
following section. The discussion of the obtained results 
and their comparison with data of Monte Carlo simula- 
tions are carried out in Sec. III. Concluding remarks are 
presented in Sec. IV. A relation between the longitudinal 
and transversal spin Green's function is checked in the 
Appendix. 

II. MAIN FORMULAS 

The Hubbard model is described by the Hamiltonian 
H = '^tii>al^aycr + ^'^ni„ni^^„, (1) 

11' (7 la 

where aj^ and aio- are the electron creation and annihi- 
lation operators, 1 labels sites of the square plane lattice, 
CT = ±1 is the spin projection, tw and U are hopping and 
on-site repulsion constants, and nio- = af^aio-. Below we 
consider the case where only the constant t for hopping 
between nearest neighbor sites is nonzero. 

In the case of strong correlations, U ^ \t\, for calculat- 
ing Green's functions it is reasonable to use the expansion 
in powers of the hopping constant. In the diagram tech- 
nique of Refs. [13 and [ij this expansion is expressed in 
terms of site cumulants of electron creation and annihi- 
lation operators. We use this technique for calculating 
the spin Green's function 

i?(lV,lr) = (T.?(r')sr(r)), (2) 

where — al^a\,^^ is the spin operator, the angular 
brackets denote the statistical averaging with the Hamil- 
tonian 

H - nig., 

1(T 

/i is the chemical potential, T is the time-ordering oper- 
ator which arranges other operators from right to left in 
ascending order of times r, and 

a\<j{T) ~ exp(7iT)aio. exp(— Hr). 

The structure elements of the used diagram technique 
are site cumulants and hopping constants which connect 
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FIG. 2; The Bethe-Salpeter equation for the sum of all four- 
leg diagrams. 

the cumulants ji^ii^ In diagrams, we denote the hopping 
constants by single directed lines. Using the diagram 
technique it can be shown that Green's function ([2]) sat- 
isfies the diagram equation plotted in Fig. [TJ In this 
diagram, after the Fourier transformation over the space 
and time variables the dual line indicates the full electron 
Green's function 

where k is the wave vector, the integer n stands for the 
fermion Matsubara frequency a;„ = {2n + 1)ttT with the 
temperature T, and (3 = T~^. The shaded circle in Fig.[T] 
is the sum of all four-leg diagrams, i.e. such diagrams in 
which starting from any leg one can reach any other leg 
moving along the hopping lines and cumulants. These 
diagrams can be separated into reducible and irreducible 
diagrams. In contrast to the latter, the reducible di- 
agrams can be divided into two disconnected parts by 
cutting two hopping lines. The sum of all four-leg di- 
agrams satisfies the Bethe-Salpeter equation shown in 
Fig. O Here the open circle indicates the sum of all irre- 
ducible four-leg diagrams. The hopping lines between the 
open and shaded circles are already renormalized here by 
the inclusion of all possible irreducible two-leg diagrams 
into these lines. These irreducible two-leg diagrams can- 
not be divided into two disconnected parts by cutting 
one hopping line.^^ As a consequence, the hopping line 
in Fig. [5] is described by the equation 

e(k,n) =ik + 4G(k,n), (3) 

where in the considered model with nearest-neighbor 
hopping we have tk = 2t[cos{kj.) + cos{ky)]. The ir- 
reducible two-leg diagrams can also be inserted in the 
external lines of the four-leg diagrams in Fig. [TJ To 
mark this renormalization we use dashed lines in that 
figure. Each of these lines introduces the multiplier 
H(k, n) — 0(k, n)/tk in the second term on the right- 
hand side of the equation in Fig. [1] Without the renor- 
malization this multiplier reduces to unity. 

As a result, the equations depicted in Figs. [1] and [2] 
read 

D{p) = -N-'TY,G{pi)G{p+Pi) 
pi 

+ Af-2y2^n(pi)n(p2)n(p + pi)H(p + p2) 

PlP2 

X ^{Pl,P+Pl,P+P2,P2), (4) 
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r{pi,p + pi,p+p2,P2) = j{pi,p + Pl,p+P2,P2) 
-N^^T^'^{pi,p + pi,p + P3,P3)Q{p3)Q{p + P3) 

P3 

X^{P3,P + P3,P + P2,P2)- (5) 

Here the combined indices p = (k, iuj^) and pj = 
(kj, iiUn^) were introduced, uji, = ivnT is the boson Mat- 
subara frequency, r(pi,p+pi,p+p2,P2) is the sum of aU 
four-leg diagrams, ^(jpi,p+pi,p+p2,P2) is its irreducible 
subset, and N is the number of sites. 

In the following consideration we simplify the general 
equations @ and ([5]) by neglecting the irreducible two- 
leg diagrams in the external and internal lines of the four- 
leg diagrams and by using the lowest-order irreducible 
four-leg diagram instead of 7(^1 , p -I- pi , p -f- P2 j P2 ) • This 
four-leg diagram is described by the second-order cumu- 
lant 

K2{t' ,t,t[,ti) = {Ta„{T')a^a{T)a_„{T[)a„{Ti))Q 



+ K,{t',t^)K^{t[,t), (6) 

where the subscript "0" of the angular bracket indicates 
that the averaging and time dependencies of the opera- 
tors are determined by the site Hamiltonian 

H\ = ^[(L//2)ni^ni _^ - /^ni^-], 

aidr) = exp(i?ir)af^exp(-iJiT), 

and the first-order cumulant 

Ki{t',t) = (Ta,(T')a,(T))o. 

All operators in the cumulants belong to the same lattice 
site. Due to the translational symmetry of the problem 
the cumulants do not depend on the site index which is 
therefore omitted in the above equations. The expression 
for K2 reads 



+e~-^°'^Ugoiini + v)goi{n2)go2{ni + n2 + v)[goi{n2 + i^) + .9oi("-i)] 
+e~^^'^Ugi2{ni + v)gi2{n2)go2{ni +n2 + v)[gi2{n2 + v) + 5i2(ni)] 
„g--Ei/3 p^^^ ^ v)gaiin2)goi{n2 + 1^) + F{n2)goi{ni + iy)goi{ni) 

+F{n2)gi2{n2 + iy)[gi2ini + iy) - goi{ni)] + F{ni + iy)gi2{ni)[gi2{n2) - goi{n2 + v)] |, (7) 

where Eq — 0, Ei — — /i, and E2 = C/— 2/i are the eigenenergies of the site Hamiltonian H\, Z = e^^"^ +26^^^^^ +e~^^^ 
is the site partition function, gij{n) = {iiUn + Ei — Ej)^^, and F{n) — 501 (") ~ 5i2("-)- 
It is worth noting that the used approximation retains the relation 

D{1't',It)=2D,{1't',It), (8) 

where 

i?,(lV,lr) = (r4(r')sf(r)) (9) 

and sf = ^ ^'^\cr^^<^ ^ component of spin. Relation ([5]) follows from the invariance of Hamiltonian ^ with 

respect to rotations of the spin quantization axisj^i The proof of Eq. ([8]) is given in the Appendix. 

Equation ([7]) can be significantly simplified for the case of principal interest U ^ T. In this case, if /i satisfies the 
condition 

e<fi<U-e, (10) 

where £ ^ T, the exponent e~^^^ is much larger than e~^^" and e^^^^. Therefore terms with e~^^° and e^^^^ can 
be omitted in Eq. ([7]) which gives 

K2ini,ni + iy,n2 + 1^,712) ^{''^(Vo + ^Sni,n2^F{ni + v)F{n2) 
-F{ni + v)gm{n2)gm{n2 + v) - F{n2)goi{ni + v)gm{ni) 

-F{n2)gi2{n2 + i^)[5i2("-i + 1^) -ffoi("-i)] - F{ni -I- 2^)512(^1) [312(^2) -5oi("2 +'^)]}- (H) 

From Eq. (O with the kernel (fTTI) it can be seen that V does not depend on momenta ki and k2. Since we 
neglected irreducible diagrams in the external lines, n(p) — 1 and in the second term on the right-hand side of 
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Eq. (|4]) the summations over ki, k2, and n2 can be carried out instantly. The resulting equation for r\^(i^,n) 
^Sn' rk('^, n + Vju' + ly, n') reads 



^'■^^{v,n) = -f^{v,n){2K2{v,n) + [a2{-v,v + n) - ai{v + n)/3(S^,o] ttkj/i(ki^) + ai{v + n)ttky2(ki^) 



where 



ai{n + i/) - a2{-iy,n + ly) + 



1 



U — iuji, 



a4{~h',72 + v) + a3{—h',n + ly) > , (13) 



(12) 



hiv, n) 



1 



1 + -F{n)F{v + n)ttk 



(14) 



, 2/i(kzy) = r^aj(i^, n)rk (1^,72), 

n 

ai{n)^F(n), a2{v,n) ^ goi{n)gai{v + n), a:i{v, n) = F{n)g i2{v + n), ai{v,n) = gi2{n) - goi{v + n). 
Multiplying Eq. (|12p by ai(z^, n) and summing over n we obtain a system of four linear algebraic equations for yi, 

Vi = h + (cj2 - Ciil35ufi)yi + Ciiy2 + Ci4j/3 + Ci3y4, (15) 

where 

h = T^ai{v,n)K2{iy,n)fk{v,n), c^j ^tt]^—'^ai{v,n)aj{-v,v + n)J^{v,n). 

n n 

System (IT5|) can easily be solved. Thus, in the used approximation the Bethe-Salpeter equation ([5]) can be solved 
exactly. In notations (fT4|) the second term on the right-hand side of Eq. (j4|) can be rewritten as 



j = ['^'^''■o(-^ "^^kyi) +ttky2] X^/kl^^, "-)ai(" + J^) + ^X!-^k(t^,")ai(n)ai(n + j/) 

- (1 - tt]^yi)^ h{f,n)a2{-iy,n + ly) + (tti^ys + jj ^ ] ^ h{'^,n)a4{-iy,n + ly) 
+ (l + ttk2/4)Xl/k(;^,n)a3(-i^,'^ + ;^)[- (16) 



In subsequent calculations we shall use the Hubbard- 
I approximatioi*^ for the electron Green's function in 
the first term on the right-hand side of Eq. In 
the used diagram technique this approximation is ob- 
tained if in the Larkin equation the sum of all irre- 
ducible two-leg diagrams is substituted by the first-order 
cumulant.^^'^"* Provided that condition pO|) is fulfilled 
the electron Green's function in the Hubbard-I approxi- 
mation reads 



III. MAGNETIC SUSCEPTIBILITY 

From the Lehmann representationi^ it can be shown 
that DzCkiy) has to be real, nonnegative, 



D,Ckiy)>0 



(18) 



G(kn) 



and symmetric, Dz(kiy) — Dz(k, —v). In view of Eq. ^ 
analogous relations are fulfilled for D{kiy). However, we 
found that condition (|18p is violated for = and some 
momentum if the temperature drops below some criti- 
cal value To which depends on the ratio \t\/U and on /i. 
As the temperature Tq is approached from above, i?(k, 0) 
tends to infinity which leads to the establishment of long- 
range spin correlations. Therefore, like in the RPA j^^i^^ 
{iujn + M)(it^n + n — U) — tk(*^ri + fJ. — U /2) wc interpret this behavior of Green's function as a tran- 

(17) sition to a long-range order. Near half-filling the highest 



iuJn + fi — U/2 
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FIG. 3: The zero- frequency magnetic susceptibility at k = 
vs. temperature at half-filling and t = — (7/4. Filled 
squares, filled and open circles are results of the Monte Carlo 
simulations,— random phase approximation, and our calcula- 
tions, respectively. 



FIG. 4: The zero-frequency magnetic susceptibility at k = 
{n,n) vs. temperature at half-filling and t = —U/A. Filled 
squares, filled and open circles are results of the Monte Carlo 
simulations,* random phase approximation, and our calcula- 
tions, respectively. 



temperature Tq occurs for the antifcrromagnetic momen- 
tum (tTjTt). Thus, near half- filling the system exhibits 
transition to the state with the long-range antifcrromag- 
netic order. 

In our calculations Tq is finite. Since we consider the 
two-dimensional model and the broken symmetry is con- 
tinuous, this result is in contradiction to the Mermin- 
Wagner theoreroi^ and shows that the used approxima- 
tion somewhat overestimates the effect of the interac- 
tion. However, it is worth noting that the value of Tq de- 
creases with decreasing the ratio \t\/U , i.e. the violation 
of the Mermin- Wagner theorem becomes less pronounced 
on enforcing the condition for which the approximation 
was developed. Notice that other approximate methods, 
including RPA'' and cluster methods, - lead also to the vi- 
olation of the Mermin- Wagner theorem. In the following 
calculations we consider only the region T > Tq. 

It was also found that for i' condition is vi- 
olated in a small area of the Brillouin zone near the T 
point. Green's function is small for such momenta and 
small negative values of D(ki/) here are a consequence 
of the used approximations. It is worth noting that the 
renormalization of internal and external hopping lines 
should improve the behavior of £'(ki/) in this region. 

To check the used approximation we shall compare our 
calculated results with data of Monte Carlo simulations^ 
on the temperature dependence of the zero- frequency sus- 
ceptibility at half-filling and on the square of the site spin 
(S^). In the usual definitions^ the susceptibility x(ki^) 
differs from D(\iv) only in a constant multiplier. For 
convenience in comparison with results of Ref. in this 
work we set 



The square of the site spin is given by the relation 



(20) 



where Eq. ([5]) is taken into account. 

The calculated zero-frequency magnetic susceptibility 
for k = and half-filling is shown in Fig. [S] Results 
obtained in Monte Carlo simulations'* and in the RPA 
are also shown here for comparison. The RPA results 
are described by the equationsi^ 



XRPA(k) 



2X0 (k) 
l-C/Xo(k)' 



(21) 



Xo(k) = -^5] 



/(ik' - m) - fih 



m) 



^k' — ik' 



X(ki.) = Dikiy). 



(19) 



where f{E) = [exp(i?/3) -f 1] *. Notice that to use 
the same scale for the susceptibility as in Ref. i our 
calculated values (fT9|) in Figs. [3] and [4] were multiplied 
by the factor 2. Also it should be mentioned that for 
T > 2\t\ = U/2 we violate condition (fTU|) : however, the 
calculated high-temperature susceptibility is in reason- 
able agreement with the Monte Carlo data. It deserves 
attention that in contrast to the RPA susceptibility which 
diverges for low temperatures the susceptibility in our 
approach tends to a finite value as it must. 

The staggered magnetic susceptibility xm is shown in 
Fig. [3) As mentioned above, in the used approximation as 
the temperature approaches Tq from above, xm tends to 
infinity which signals the establishment of the long-range 
antifcrromagnetic order. For parameters of Fig. [3] Tq « 
0.64|t|. The transition temperature Tq is finite; however, 
for the considered range of parameters 4|t| < U < 16|t| 
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FIG. 5: The square of the site spin (S^) vs. tempera- 
ture at half-filling. Filled symbols are data of Monte Carlo 
simulations,^ open symbols are our results. Squares and cir- 
cles correspond to the cases t = ~U/8 and t = — f//4, respec- 
tively. 



it is always lower than the respective temperature in the 
RPA. Accordingly our calculated values of xm in Fig. 0] 
are closer to the Monte Carlo data than the RPA results. 

The temperature variation of the square of the site 
spin, Eq. (jSO]) , is shown in Fig. \5\ together with the data 
of Monte Carlo simulations.* As might be expected, the 
results for the smaller ratio \t\/U more closely reproduce 
the data of numerical simulations. For t ~ — C//8 our 
calculations replicate the Monte Carlo data for T > \t\ 
and the difference between the two series of results is 
less than 10 percent. This difference is at least partly 
connected with the simplification made above when irre- 
ducible two-leg diagrams were dropped from internal and 
external lines of the four-leg diagrams. The difference be- 
comes even smaller if in accord with the Mermin- Wagner 
theorem Tq is set as the zero of the temperature scale and 
our calculated curve is offset by this temperature to the 
left. On approaching Tq our approximation becomes in- 
applicable for calculating (S^) - it starts to grow rapidly 
and exceeds the maximum value |. 

The concentration dependence of (S^) near half- filling 
is shown in Fig. [51 The range of the electron concentra- 
tion n = '^g.{nia) which corresponds to the chemical po- 
tential shown in this figure spans approximately 0.8 — 1.2 
for t — —U/8. As would be expected, (S^) decreases 
rapidly with the departure from half-filling. 

The momentum dependence of the zero-frequency sus- 
ceptibility at half-filling and its variation with tempera- 
ture are shown in Fig. [71 At half-filling the susceptibility 
is peaked at the antiferromagnetic wave vector (7r,7r). 
For temperatures which are only slightly higher than Tq 
the peak intensity is large [Fig. [7l (a)] which leads to a 
slow decrease of spin correlations with distance and long 
correlation lengths (see below) . With increasing temper- 
ature the peak intensity of the susceptibility decreases 



' 2 ' 4 ' 6 ' 8 

FIG. 6: The square of the site spin (S^) vs. the chemical 
potential for t = -17/8 and T = \t\. 



rapidly [Fig. [71 (b) and (c)] which results in a substantial 
reduction of the correlation length. In this case for dis- 
tances of several lattice periods the spin correlations are 
small, nevertheless they remain antiferromagnetic. 

The situation is changed with the departure from half- 
filling. The zero-frequency susceptibility for different 
electron concentrations is shown in Fig. [SI The values 
of the concentration which correspond to parts (a) to (c) 
are n « 0.94, 0.88, and 0.81, respectively. Notice the 
rapid decrease of the peak intensity of the susceptibility 
with doping [cf. parts (a) in this and the previous fig- 
ure]. Starting from n « 0.94 the susceptibility becomes 
incommensurate - the maximum value of the susceptibil- 
ity is not located at (tt, tt) - and the incommensurability 
parameter, i.e. the distance between (tt, tt) and the wave 
vector of the susceptibility maximum, grows with depar- 
ture from half-filling. It is interesting to notice that for 
n < 1 the zero-frequency susceptibility diverges when the 
temperature approaches some critical temperature in the 
same manner as it does at half- filling. For t = —U/8 and 
0.94 < n < 1 the divergence first occurs at (7r,7r), while 
for smaller electron concentrations it appears at incom- 
mensurate wave vectors. For n < 1 the value of the 
critical temperature is less than Tq - the temperature 
at which the transition to the long-range order occurs 
at half-filling. The critical temperature decreases with 
decreasing n. If in accord with the Mermin- Wagner the- 
orem we identify Tq with zero temperature we have to 
conclude that for n < 1 the system undergoes a virtual 
transition at negative temperatures, while for T > it is 
governed by short-range order. In view of the particle- 
hole symmetry analogous conclusions can be made for 
n > 1. 

Analyzing equations of the previous section it can 
be seen that the momentum dependence of the zero- 
frequency susceptibility is mainly determined by the mul- 
tiplier 2/1 (k, v — 0) in the first term on the right-hand side 
of Eq. (fTB]) . At half-filling the susceptibility is commen- 
surate, since this term is peaked at (tt, tt) and diverges at 
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FIG. 7: (Color online) The zero-frequency magnetic suscepti- piG. 8: (Color online) The zero-frequency magnetic suscepti- 
bility at half-filling for f = -17/8 in a quadrant of the Brillouin bility for t = -U/8 and T = 0.06C/ in a quadrant of the Bril- 
zone. (a) T = 0.06(7, (b) T = 0.1(7, and (c) T = 0.2(7. louin zone, (a) /i = 0.2(7, (b) /i = 0.15(7, and (c) fj. = 0.1(7. 



this momentum when T ^ +7o, as the determinant of 
the system P3|) vanishes. At departure from half-filling 
the behavior of yi is governed by the term bi in this 
system. The term contains the sum 



r^a?(0,n)/k(0,n) 



1 



1 



(22) 

1. For half- 



where F{n) = — C/[(«a;„ -I- ^){iujn + /i — (7)] 
filling the sum has a maximum at (tt, tt), however with de- 
parture from half-filling the maximum shifts from (tt, 7r~ 
and the susceptibility becomes incommensurate. 



Together with the zero-frequency susceptibility the 
imaginary part of the real-frequency susceptibility, 



x"(kijj) = ImZ?(k, a; + 177), 77 



-0, 



(23) 



becomes also incommensurate. This quantity is of spe- 
cial interest, because it determines the dynamic structure 
factor measured in neutron scattering experiments 1^ To 
carry out the analytic continuation of D{kv) to the real 
frequency axis an algorithm.^^ based on the use of Pade 
approximants can be applied. In this calculation 300 val- 
ues of Z?(ki^) at equally spaced imaginary frequencies in 
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FIG. 9: (a) The momentum dependence of x"(ka)) along the 
edge [sohd hne, k = (tt, k)] and diagonal [dashed line, k — 
(k, k)] of the Brillouin zone for t = -0.11(7, lu = 0.002(7 and 
n ~ 0.88). (b) The momentum dependence of x"(k(x;) along 
the zone edge for n ~ 0.88 (solid line), n ~ 0.94 (dashed line), 
and n = 1 (dash-dotted line), t = -0.11(7 and lo = 0.002(7. 
(c) The dispersion of maxima in x"(ktj) along the zone edge 
for t = -0.11(7 and n 0.88. 



the upper half-plane were used. The obtained dependen- 
cies of the susceptibility on the momentum for a fixed 
transfer frequency uj and the dispersion of low-frequency 
maxima in x" are shown in Fig. [51 The susceptibility is 
shown in the first Brillouin zone and can be extended to 



the second zone by reflection with respect to the right 
y axis. As seen from Figs. M (a) and (b), with depar- 
ture from half-filling x"(kti;) becomes incommensurate 
and the incommensurability parameter grows with in- 
creasing 1 — n. 

This behavior of the susceptibility x"(kw) in the Hub- 
bard model resembles the low-frequency incommensurate 
magnetic response observed by inelastic neutron scatter- 
ing in lanthanum cupratesii^ In these crystals, the in- 
commensurability is observed both in the normal and 
superconducting states. For small transfer frequencies lo 
the maxima of the susceptibility are located on the edge 
of the Brillouin zone. For the parameters of Fig. [5] (a) 
our calculated susceptibility is also peaked on the zone 
edge. However, for other parameters the susceptibility 
on the diagonal may be comparable to that on the zone 
edge. This uncertainty in the position of the suscepti- 
bility maxima may be connected with errors introduced 
in the calculation results by the procedure of analytic 
continuation to real frequencies. 

In experiment, for small lo the incommensurability pa- 
rameter grows with the hole concentration 1 — n in the 
range 0.04 ^ 1 — n < 0.12 and saturates for its larger val- 
ues. This behavior of the incommensurability parameter 
is reproduced in our calculations [see Fig. ^ (b)] and its 
values are close to those observed experimentally. For a 
fixed hole concentration the incommensurability param- 
eter decreases with increasing to and at some frequency 
LOr the incommensurability disappears and the suscepti- 
bility x"(kw) appears to be peaked at the antiferromag- 
netic momentumi^ The same behavior is observed in the 
Hubbard model [see Fig. [n](c)]. In lanthanum cuprates 
for the hole concentrations 1 — n « 0.12 the frequency 
u!r « 50 meV. In Fig. [9] (c) we chose parameters so that 
LOr was close to this value (for the superexchange con- 
stant J = W^/U w 0.15 eV and t = -0.11(7 we find 
U = 3.1 eV,t = 0.34 eV, and lo^ = 44 meV). Notice that 
like in experiment LOr decreases with decreasing 1 — h. 

A similar incommensurability is observed in 
YBa2Cu307_yj^ however, in this case due to a 
larger superconducting temperature and gap the mag- 
netic incommensurability is usually observed in the 
superconducting state and the low-frequency part of 
the susceptibility is suppressed. As follows from the 
above discussion, in the Hubbard model the magnetic 
incommensurability is a property of strong electron cor- 
relations. The similarity of the mentioned experimental 
and calculated results gives ground to consider these 
strong correlations as a possible mechanism of the low- 
frequency incommensurability observed in experiment. 
A similar mechanism was observed for the related t-J 
model in Ref. [13. 

In experiment (^i^^ for frequencies oj > lo^. the suscepti- 
bility x"(\iLo) becomes again incommensurate such that 
the dispersion of maxima in x"(kci;) resembles a sand- 
glass. The most frequently used interpretations of this 
dispersion are based on the picture of itinerant electrons 
with the susceptibility calculated in the RPA^^ and on 
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FIG. 10: The susceptibility x"(kt^) for k = (vr, n), t = -17/8, 
T = 0.06[/, and = 0.2[/'(n f» 0.94). 

the stripe picture/^i^^ In Ref . 12^ the sandglass dispersion 
was obtained in the t-J model in the regime of strong 
electron correlations without the supposition of the ex- 
istence of stripes. In this latter work the part of the 
sandglass dispersion for uj > uJr was related to the dis- 
persion of excitations of localized spins. Similar notion 
was earlier suggested in Ref. [13. In our present calcula- 
tions we did not obtain this upper part of the dispersion, 
since the used approximation does not describe the ap- 
pearance of localized spins. A typical example of the 
frequency dependence of the susceptibility x"(kw) which 
up to the multiplier tt"^ coincides with the spin spec- 
tral function is shown in Fig. [TOl The susceptibility usu- 
ally contains several maxima one of which is located at 
u! <^ U, while others are placed at frequencies of the 
order of U. Since the localized spin excitations have fre- 
quencies in the range < < 2 J where J = /U <^ U, 
the former maximum could be taken as a signal for such 
excitation. However, the intensity of the maximum usu- 
ally grows with temperature and with departure from 
half-filling. This indicates that the maximum is more 
likely due to a bound electron-hole state in which both 
components belong to the same Hubbard subband, while 
in the high-frequency maxima the components belong to 
different subbands. 

In connection with the Nagaoka theoremSS it is of inter- 
est to investigate the tendency towards the establishment 
of the ferromagnetic ordering with departure from half- 
filling. For a finite U this problem was investigated by 
different analytical methods^ '^^i^^'^" and by Monte Carlo 
simulations.^ Our results for the spin-spin correlator, 

^4-^o> = §Ecos(kL)i5('''^)' (24) 

as a function of the distance between spins are shown 
in Fig. [TT]for different parameters. Figure[Tl] (a) demon- 
strates the short-range antiferromagnetic order at half- 
filling for a temperature which is slightly above Tq (as dis- 




FIG. 11: The spin-spin correlator (s^Sq) for L = {Lx,0) and 
t = -17/8. (a) r = OmU, fi = 0.5(7, (b) T = Q.125f/, 
fi = 0.5U, and (c) T = 0.06(7, fi = 0.1(7 (n « 0.81). Insets in 
(b) and (c) demonstrate the same data as in the main plots 
in a larger scale. 

cussed above in connection with Fig. [5l for such temper- 
atures the value of (sq Sq ) is somewhat overestimated by 
the used approximation). Figure [TT] fb) corresponds also 
to half-filling to somewhat higher temperature. In this 
case the correlations are still antiferromagnetic though 
they are characterized by a correlation length which is 
much shorter than that in Fig. [11] (a). We have found 
that the correlation length diverges when T — > Tq which 
indicates the transition to the long-range antiferromag- 
netic order. Similar weak antiferromagnetic correlations 
were also obtained for moderate departures from half- 
filling. Figure [Tl] (c) corresponds to the lowest filling 
h « 0.81 which is allowed by condition PH)) for the given 
ratio J7/|t|. According to the mean-field theory"* and the 
generalized RPAiS in this case the system has a ferro- 
magnetic ground state. As seen from Fig. [Tl] (c), our 
calculated spin-spin correlations are still antiferromag- 
netic even for nearest neighbor spins. This result is in 
agreement with Monte Carlo simulations'* carried out for 
the same parameters. Analogous result was also obtained 
for U = 16\t\. However, a tendency for the establishment 
of ferromagnetic correlations can also be seen from the 
comparison of Figs. [11] (a) and (c) - the antiferromag- 
netic spin correlation on nearest neighbor sites becomes 
smaller with doping. For larger ratios of U/\t\ we can 
ascertain that the correlation changes sign and becomes 
ferromagnetic. In particular, it happens at U/\t\ = 25 
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and n w 0.77. For these parameters condition (fTO|) is 
still fulfilled. 



IV. CONCLUDING REMARKS 

In this article we investigated the magnetic suscepti- 
bility of the two-dimensional repulsive Hubbard model 
using the diagram technique developed for the case of 
strong electron correlations. In this technique the power 
series in the hopping constant is used. At half-filling the 
calculated temperature dependence of the zero-frequency 
susceptibility reproduces adequately key features of re- 
sults of Monte Carlo simulations. The uniform suscep- 
tibility tends to a finite value for vanishing tempera- 
ture. The staggered susceptibility diverges with decreas- 
ing temperature which signals the establishment of the 
long-range antiferromagnetic order. The transition tem- 
perature is finite which indicates the violation of the 
Mermin- Wagner theorem. However, the transition tem- 
perature is always lower than the analogous tempera- 
ture in the RPA. Besides, the transition temperature de- 
creases with the decrease of the ratio \t\/U of the hopping 
constant and the on-site repulsion, i.e. the violation of 
the Mermin- Wagner theorem becomes less pronounced 
on enforcing the condition for which the approximation 
was developed. For small ratios \t\/U the calculated 
square of the site spin differs by less than 10 percent from 



the data of Monte Carlo simulations. Also in agreement 
with Monte Carlo results we found no evidence of ferro- 
magnetic correlations in the considered range of electron 
concentrations 0.8 ^ n ^ 1.2 for the repulsion param- 
eters 8|i| < U < 16\t\. However, for larger U/\t\ and 
1 — n| « 0.2 the nearest neighbor correlations become 
ferromagnetic. In the case U = 8\t\ for 0.94 ^ n ^ 1.06 
the zero-frequency susceptibility and the imaginary part 
of the susceptibility for low real frequencies are peaked 
at the antiferromagnetic wave vector (tTjTt). For smaller 
and larger concentrations these susceptibilities become 
incommensurate - momenta of their maxima are shifted 
from (tt, tt) - and the incommensurability parameter, i.e. 
the distance between (tt, tt) and the momentum of the 
maximum susceptibility, grows with departure from half- 
filling. With increasing the transfer frequency the incom- 
mensurability parameter decreases and finally vanishes. 
This behavior of the susceptibility in the strongly corre- 
lated system can explain the observed low-frequency in- 
commensurate response in the normal state of lanthanum 
cuprates. 



Acknowledgments 

This work was partially supported by the ETF grant 
No. 6918 and by the DFG. 



APPENDIX 



In this Appendix we prove the symmetry relation ([5]) . In the zeroth order of the perturbation expansion for Green's 
function ^ we find 



Di°^l'T',lT) = -^a(T'\K'2{T'cj',T'a',Ta,Ta)5iv + Ki{t' t')Ki{tt) - Ki{t' T)Ki{TT')5iv5aa' 
= ^6u' ^K'2{T'a,T'a,Ta,Ta) - /^alr'o-; r'cr; r, -ct; r, -a) - Ki{t't)Ki{tt') 



(A.l) 



where we took into account that the first-order cumulant Ki{t't) does not depend on a and therefore the second 
term in the sum vanishes. Up to the multiplier i the last term on the right-hand side of Eq. (jA.ip coincides with 
the respective term in the expansion for Green's function In the used diagram technique Ki describes the bare 
electron Green's function. Therefore that term contributes to the electron bubble shown in Fig. [T] From the higher- 
order terms it can be seen that inclusion of irreducible two-leg diagrams into the two bare Green's functions of that 
term retains the one-to-one correspondence between terms of the bubble diagrams in D and Dz and the additional 
multiplier ^ in the terms of D^- 

The second-order cumulant K2 in Eq. (jA.l[) is defined as 

K'2iT'cT,Ta,T[cTi,Ti<7i) = (Ta, (t' (t (rOa.i ) ) - (t' , r)if 1 (t{ , Ti ) -f i^i (t' , Ti )Ki (t{ , r)5,,, . (A.2) 

This definition is more general than Eq. ^ - the latter is obtained from Eq. (|A.2p if we set tri = —a and interchange 
annihilation operators in K2. After the Fourier transformation we find 



K2{nia;ni + iy,cr;n2 + iy,ai;n2<7i) = Z ^|/3 {S^oSaai ~ Snini)^ 



F{ni + v)F{n2) 



11 



-(^o-.-criG ^"^Ugoiini + i^)goi ("2)302 (n-i + "2 + i^) [.goi("-2 + v) + gm{ni)\ 
-5a-aie^^^^Ugi2{ni + ^^)gl2(r^2)go2(?^l + ^2 + i^) [.912(^2 + i^) + gl2(?^l)] 
+5a-aie^^^^^ F{ni + J^)5oi("-2)ffoi("2 + i') + F{n2)goi{ni + i').9oi("i) 

+F(n2)5i2(ri2 + i^) [.9i2(n.i + i^) - goi('^i)] + -^("-1 + i^)5i2(ni) [gi2(n2) - .9oi("2 + t^)] |, (A. 3) 

where the notations are the same as in Eq. ([7]). From these two equations it can be seen that 

K2{ni,ni + z/, ^2 + I', 712) = K2{nia; ni + v, cr; n2 + cr; ri2cr) — K'2{ni<j] ni + cr; ri2 + — cr; ^2, -cr) (A. 4) 

and the analogous equation is fulfilled for the Fourier-transformed quantities. Thus, zeroth-order terms in the expan- 
sions for D and coincide up to the factor i. 

The next terms in the considered expansions for D and Dz contain two second-order cumulants and appear in the 
second order. These terms read 



i?(2)(lV',lr) = - 11^^ dndT2tn'tnK2{T',T',T,,T2)K2{T2,n, 



T, r), 



Di^)(l'rMr) 



dTidT2tu'triK2{T'(T' , t'ct', Ti(Ji,T2ai)K2{T2ai,Tiai,Ta, Ta) 



(A.5) 
(A.6) 



Using twice relation \KA\ in Eq. (|A.6p one can see that i:>(2) ^ 2Df\ Analo gous equations for higher order terms 
can be proved in the same manner. Thus, relation ([8]) is fulfilled. 
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